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Abstract 

The construction of measurements suitable for discriminating signal components pro- 
duced by phenomena of different nature is considered. The required measurements should 
be capable of cancelling out those signal components which are to be ignored when fo- 
cussing on a phenomenon of interest. Under the hypothesis that the subspaces hosting 
the signal components produced by each phenomenon are complementary, their discrimi- 
nation is accomplished by measurements giving rise to the appropriate oblique projector 
operator. The subspace onto which the operator should project is selected by non-linear 
techniques in line with adaptive pursuit strategies. 

PACS numbers: 02.30.Mv, 02.60.Gf, 02.30.Sa, 83.85.Ns, 95.75.Tv, 95.75.Pq, 95.75.Fg 

1 Introduction 

The word signal is frequently used to refer to a physical carrier convening information about 
some phenomenon. We adopt such terminology and further refer to the process of transforming 
a signal into a number (within the corresponding units) as a measurement. An appropriate 
mathematical setting for this description is to consider that a signal is an element of some 
vector space and a measurement a functional transforming the vector into a scalar. In this 
effort we discuss the design of measurements in relation to the following problem: Assume that 
a signal /, represented as an element of an inner product space H, arises by the superposition 
of two components, f\ and f 2 , each of which is produced by a particular phenomenon and such 
that f% G <Si and / 2 G S 2 , where «Si and S 2 are disjoint subspaces of TC, i.e. <Si fliS>2 = {0}. This 
condition implies that the superposition / = fx + f 2 is unique. The matter to be addressed 
here concerns the construction of the appropriate measurements allowing us to discriminate 
the component, say fx, from the available signal / and the knowledge of Sx and S 2 . Under 
the condition Sx PI S 2 = {0} the problem has a straightforward 'theoretical' solution, since 
the component fx can be extracted from / by an oblique projection onto Sx and along S 2 
[1,2]. Unfortunately, even when theoretically the condition Sx fl S 2 = {0} is satisfied, if 
the subspaces Sx and S 2 are not well separated, the construction of the corresponding oblique 
projector becomes ill posed. Consequently, the signal splitting can not be achieved by numerical 
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calculations in finite precision arithmetics. This is the situation we are concerned with. We 
assume that the given subspaces Si and S% are 'theoretically' disjoint, but close enough to yield 
an ill posed problem. 

Our proposal for the numerical realization of the phenomenon discrimination is focussed on 
the search of a subspace of the given Si, where a class of signals is considered to lie. It will 
be assumed throughout the paper that the class of signals to be considered is K-sparse in a 
spanning set for <Si. By this we mean that given a spanning set for Si, the corresponding linear 
superposition of a signal has at most K nonzero coefficients. The i^-value should be less than 
or equal to the dimension of the subspace of Si such that the construction of measurement 
vectors giving rise to an oblique projection onto itself, and along <S 2 , is well conditioned. This 
assumption is quite realistic, considering that in practice there is often a lack of complete 
knowledge on the actual subspace Si and to be on the safe side one may overestimate it. 
However, the assumption does not make the problem much easier to solve. Indeed, the problem 
of subspace selection is in general a combinatorial problem, whereby an exhaustive search of 
possibilities is in general intractable. The approach we propose in this Communication evolves 
by step wise optimal selection and is in line with the adaptive greedy approximation termed 
Matching Pursuit. Such a technique, which appeared first in the statistic literature [3], has 
been extended in the area of signal processing to several greedy strategies [4-9] being currently 
of assistance to a range of disciplines, including physics [10-12]. In particular, we revise and 
extend the Oblique Matching Pursuit (OBMP) approach which has been recently proposed in 
relation the above described problem of signal discrimination [13]. 

The paper is organized as follows: In Section [2] we introduce the mathematical setting for 
signal representation to be adopted here and in Section [3] with discuss the construction of 
oblique projectors. Section H] highlights the importance of the search for sparse representations 
in the construction of oblique projectors for phenomenon discrimination. The proposed strategy 
is discussed in Section \5\ and illustrated in Section E] by two numerical simulations: l)The 
cancellation of impulsive noise from the register of a system of harmonic oscillators and 2) the 
separation of a spectrum from blackbody radiation background. The conclusions are presented 
in Section [3 

2 Mathematical framework 

Regardless of the informational content of a signal, we deal with it mathematically by consider- 
ing it as an element of an inner product space 7i. Thus, adopting Dirac's notation, we represent 
a signal / as a ket |/) and the corresponding dual as a bar (/|. Accordingly, the square norm 
|||/)|| 2 is induced by the inner product that we indicate as (f\f)- For our present purpose we 
further assume that all the signals of interest belong to some finite dimensional subspace V G 7i 
spanned by a finite set {\vi) G 7Y}*£ 1 . Consequently, for every signal |/v) G V there exists a set 
of numbers {q}*^ which allow us to express the signal as the linear superposition 



In the jargon of signal processing the above expansion is called atomic decomposition and the 
vectors in the decomposition are called 'atoms' [14]. In applications where an economical signal 
representation is important, the goal is to construct decompositions involving as few terms as 
possible. For this end the atoms are selected from a large and, in general redundant, set called a 
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dictionary. If the number of M-terms in ([I]) is satisfactorily small in relation with the particular 
application, the decomposition is said to be sparse. 

Even when we think of a signal as an abstract object in an inner product space, for process- 
ing tasks we need a numerical representation of such an object. The process of transforming 
the signal into a number is refereed to as a measurement or sampling. Since we restrict consid- 
erations to linear measurements, we represent them by linear Junctionals. Thus, making use of 
Riesz' theorem [15] we can express a linear measurement as m = (w\f) for some \w) e 7i. Hence, 
considering M measurements rrii, i — 1, . . . , M, each of which is obtained by a measurement 
vector \wi) we have a numerical representation of the ket |/) as given by 

m i = (w i \f), i = l,...,M. (2) 

The representation of measures as in ([2]), which have been used in physics for many years, 
has started to become popular within other disciplines through the theory of Compressed 
Sensing [16-20]. 

The problem of reconstructing the whole information content of a signal from a numerical 
representation has been extensively studied for the last thirty years from a number of different 
points of view. The diversity of available approaches is specially helpful when the reconstruction 
is to be achieved on the basis of incomplete information. In particular, the above mentioned 
theory of Compressed Sensing has produced strong theoretical results with regard to the recov- 
ery of a signal, assumed to be sparse in some orthonormal basis, from a number of non adaptive 
measurements which can be significantly less than the dimension of the signal subspace [16-20]. 
Here we focus on the particular problem involving the reconstruction of a signal from a set of 
linear measurements as given in (T5]), with no further numerical calculations other than using 
these numbers in the expansion (CQ). In other words, we wish to use the measurements as coeffi- 
cients in the linear combination The question then arises as to which are the conditions to 
be requested for the measurement vectors \wi) to produce the corresponding numbers allowing 
the reconstruction of the signal |/) as 

M 

l/v> = £>>H/>. (3) 

i=l 

A major consequence of working under the assumption that the signal of interest, \fv), be- 
longs to a finite dimensional subspace, V, is the lack of uniqueness of the measurement vectors 
{\ w i)}t£i> even when the spanning set {{vi)}^ is linearly independent. This statement appears 
clearly from the following observation. 

Let us denote Y^iLi \ v i)( w t\ as an operator, E, so as to recast ([3]) in the fashion 

Lfv> = %>. (4) 

This equation tells us that the measurement vectors {{wi)}^ should be such that operator E is 
a projector onto V. Indeed, the operator E is a projector if and only if it is idempotent [15,21], 
i.e., E 2 = E. Consequently, as discussed below, the projection is onto the range of the operator, 
7Z(E), and along its null space M{E). 

Denoting by T> the domain of E we recall that 

n{E) = {|/), such that |/) = E\g), for some|#) G V}. 
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Thus, for E an idempotent operator and for |/) G 71(E), we have E\f) = E 2 \g) = E\g) = \f). 
This implies that E behaves like the identity operator for all |/) G 71(E), regardless of M(E), 
which is defined as 

M(E) = {\g), such that E\g) = 0, \g) G V}. 

It is now clear that to reconstruct a signal |/) G V by means of ([3]) the measurement vectors 
{\ w i)}iL\ should give rise to an operator ^i=i \ v i)( w i\i which must be a projector onto V. It 
is appropriate to point out that the required operator is not unique (even if the spanning set 
is linearly independent) because there exist many projectors onto V having different 
J\f(E). Consequently, for reconstructing signals in the range of the projector its null space can 
be chosen arbitrarily. Nevertheless, the null space, and therefore the particular measurement 
vectors, become crucial when the projector is to be applied on signals outside its range. It 
follows then that the measurement vectors {|wi)}i=i can be tailored for a particular purpose. 
Such a degree of freedom will be indicated hereforth by using two subscripts for representing 
a projector. We adopt the notation E vw ± to indicate a projector onto the subspace V and 
along the subspace W x . The particular case E vv ±, where V -1 is orthogonal to V, corresponds 
to an orthogonal projector and we use the special notation Py to denote such a projector. The 
orthogonal projector is popular in approximation techniques because if a signal |/) is to be 
approximated by a signal |/y) G V the choice |/y) = Pv\f) is known to yield the unique signal 
in V minimizing the distance |||/) — |/y)||. However, as discussed below, if one is interested 
in discriminating signal components produced by phenomena of different nature an alternative 
selection of the subspace W -1 is required. When W -1 is not orthogonal to V the projector E vw ± 
is referred to as an oblique projector. 

Let as assume for instance that a signal |/) is the superposition of two signals, |/) = |/y) + 
I/w-l), each component being produced by a different phenomenon we wish to discriminate. Let 
us assume further that we can model the subspaces V and W ± hosting each signal component 
and such subspaces are disjoint, i.e. V fl W -1 = {0}. Thus we can obtain \f v ) from |/), by an 
oblique projector onto V and along W L . The projector will map to zero the component |/w±) 
to produce 

\f v ) = E vw x\f). 

In the next section we discuss the construction of measurement vectors giving rise to 

the desired projector. 

3 Constructing measurement vectors for discrimination 
of signal components 

Given two disjoint subspaces V and W" 1 , in order to provide a prescription for constructing the 
projector E vw ± one can proceed as follows. Firstly we define S as the direct sum of V and 
W" 1 , which we express as 

s = v®w ± . 

Let W = (W -1 )" 1 " be the orthogonal complement of W -1 in S. Thus we have S = V © W -1 = 
W © W ± , where the operation ©- 1 indicates the orthogonal sum, which refers to the direct 
sum of orthogonal subspaces. 

Considering that is a spanning set for V, a spanning set for W is obtained as 

K) = k) - P\vA v i) = Avk), i = l,...,M. (5) 
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Denoting as the standard orthonormal basis for C M , we define the operators V : C M 

V and U : C A/ -> W as 

M M 

i=i i=i 

Thus the adjoint operators U* and V* are expressed as 

M M 



v* = Y,m< ft* =!>><«« 



i=l i=l 

Notice that Pyv^ = # and f>*P w = f>* hence, G : C M -> C M defined as: 

G = U*V = U*U 

is a self-adjoint operator. Its matrix representation being given by the elements (i\G\j) = 
(ui\vj) = (ui\uj), i,j = 1,...,M. 

Remark 1. It is appropriate to stress that 

• Operators V and U are given in terms of spanning sets for the spaces V and W, respec- 
tively, and any such spanning set can be used. 

• The condition V D W -1 = {0} implies that the dimension ofV is equal to the dimension 
ofW. Hence, provided that the spanning set {{vi}}^ is linearly independent, operator 
G has an inverse. Nevertheless, the independence 0/{|%}}^i is not a requirement and 
therefore an inverse for G need not exist. For the sake of generality we shall use G\ 
which indicates a pseudo-inverse of G. 

The oblique projector operator onto V and along W is given as [22] 

E vw ± = V&U*, (6) 



or, equivalently, as 



with 



M 



Eyw^- =^2\vi)(wi\, (7) 



i=i 

M 

\w t ) = U&\t)=J2 = \^)(j\&W. 



It is actually straightforward to verify that E vw ± given in (JTj) satisfies the required properties. 
Namely, ££ w± = £ vw ±, E vw ±\f v ) = \fv) for all |/ v ) G V, and E vw x\g) = or all \g) G W ± . 

Remark 2. The construction of an oblique projector is similar to that of an orthogonal one. The 
difference being that in general the subspaces span^Vj)}^ = V and span^Wi)}^ = W are 
different. For the special case {|fi)}|£i = {|wj)}|£i we have span^Wi)}^ = span^Vi)}^ = V , 
thereby the projector is self adjoint and, consequently, an orthogonal projector onto V along V ± . 
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We are already in a position to extract the component |/v) from |/) by the simple operation 
\fv) — VGW*\f). However, the correct discrimination of the signal components is successful 
provided that the subspaces V and W 1 - are well separated. Unfortunately, this is not always the 
case and the construction of the necessary projector may generate an ill posed problem. In spite 
of the fact that 'theoretically' V H W ± = {0}, numerical errors, due to the existence of small 
eigenvalues values of the operator G, may cause the failure to find the unique signal splitting that 
theoretically one should expect. Nevertheless, the correct separation is still possible, provided 
that the signal |/y) admits a sparse representation in some spanning set for V. In order words, 
one could succeed in extracting \fv), provided that it is well represented in a subspace Vk C V 
inducing a subspace Wk C W (satisfying Vk + VV 1 - = WV ©- 1 W -1 ) where the computation of 
the measurement vectors is well posed. If this is the case, the problem of designing measurement 
vectors for discriminating signal components can be addressed as the problem of finding the 
subspace Vk where \f v ) is well represented. Unfortunately the search for the subspace Vk is 
in general intractable. Indeed, let us assume that is an spanning set for V and Vk is 

spanned by K elements of such a set. Even possessing this knowledge, the problem of finding 
the right subspace by exhaustive search would be a combinatorial problem: Out of a set of 
cardinality M there exist possible subsets of cardinality K. As already mentioned, we 
shall not look for the sparsest representation but make the search of the appropriate subspace 
tractable by means of recursive greedy pursuit strategies, which are only step wise optimal. 
Before discussing our approach some considerations are in order. 

4 Getting ready for the search 

In this section we highlight some properties that will be of assistance in the next section, where 
we will present our strategy for the search of the sparse representation achieving the desired 
signal discrimination. The goal is to avoid the computation of the measurement vectors in the 
whole subspace. Instead, we strive to find the subspace Vk C V, where the signal component 
one wants to extract from a signal |/) is assumed to lie. We work under the hypothesis that 
the subspace W 1 - is given and fixed. Furthermore, V D W -1 = {0}, which implies that there 
exists a unique solution for the signal splitting. The problem we need to address arises from the 
fact that, if the subspaces V and W ± are not well separated, the numerical calculation of the 
measurement vectors is not accurate (due to the numerical operations being carried out in finite 
precision arithmetic). As a consequence, the representation of the corresponding projector fails 
to produce the correct signals separation. This effect is very much magnified if the data are 
affected by errors no matter how insignificant those errors are. 

Assuming that we are able to accurately compute in finite precision arithmetic r measure- 
ment vectors, we could attempt to single out a signal belonging to a subspace spanned by at 
most r vectors (i.e. we could attempt to separate from |/) a signal expressible as in ([3]) but at 
most with r nonzero coefficients). However, as discussed above, even possessing this knowledge 
about the sought signal the problem of finding the right subspace by exhaustive search is not 
affordable. Hence, an adaptive greedy strategy for the subspace selection, given a signal, was 
advanced in [13]. Before revising and extending that strategy we need to recall two relevant 
properties of oblique projectors. 

Property 1. The oblique projector E vw ± satisfies Py\>E vyv i_ = Pyv- 

Proof. It readily follows by applying P w on both sides of (jBD or (J7J). Since AvK) — \ u i) an d 
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(ui,Vj) = (ui,Uj), one has 

M 

P w E vw x = \m)(v>i\ = u&u* = Av- (9) 

□ 

Property 2. Given a signal \f) inV + W 1 - = W © W -1 , the only vector \g) G V satisfying 

Pw\f) = Pw\g) (io) 

m l#) = E vw ±\f). 

Proof. If |g) = E vw ±\f) (TTDT) trivially follows from Property [TJ Let us assume now that there 
exists \g) G V such that (PUD holds. Then Av(|/> - = 0, i.e., (|/) - \g)) G W^. Hence 
^'vvv- L (l/) ~~ Ifl 1 )) = an d, since \g) G V, this implies that E vw ±\f) = \g). □ 

Let us suppose that Vk = span{|vj)}f =1 is given and the spanning set is linearly independent. 
Assuming that V^nW 1 = {0} we guarantee that the set of vectors {\ui)}^ =1 , with \ui) given in 
([5]), is also linearly independent. Consequently the dimension of Vk is equal to the dimension 
of Wfc = span{|wj)}f =1 = span{|u>f)}f =1 . We use now a superscript A; to indicate that the 
measurement vectors {|u>f)}f =1 span Wk- Hence these vectors give rise to the oblique projection 
of a signal |/), onto Vk and along W" 1 , as given by: 

k k 

E Vk wAf) = X>x«>?i/> = E c ^>- ( n ) 

i=l i=l 

It is clear from (|TT|) that if the atoms in the atomic decomposition were to be changed (or some 
atoms were added to or deleted from the decomposition) the measurement vectors {|wf)}f =1 , 
and consequently the coefficients {cf}f =1 in (fTTj) . would need to be modified. The recursive 
equations below provide an effective way of implementing the task. 



Forward/backward adapting of measurement vectors 

Starting with \w\) = ip^up ? an d \ux) as in (j3j), the measurement vectors {\w^ +1 }}^ can 
be recursively constructed from {|u>f)}^ =1 as follows [2]: 

K fc+1 ) = k*>-l«4ti><«w-ik?>, i = i,...,k (12) 

= ii } lk+1 }\\2 i l7fc+i> = l«fc+i) - Avn|wfc+i)> ( 13 ) 
Ill7fc+i/ll 

where Py\> k is the orthogonal projector onto Wk = span{|Mj)}^ =1 . We note that, since = 
Pyv\vk+i) and Pwl^f) = \ w f), (O can also be written as 

k fc+1 > = H) - i = i,...,k. (i4) 

It follows from the above equations that when incorporating a linearly independent atom |ffc+i) 
in the atomic decomposition (TTTT) . the coefficients can be conveniently modified according to 
the recursive equations 

<&l = «il/>, (is) 

= K fc+1 |/) = cf-c^K fc K +1 ), i = l,...,k. (16) 
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Conversely, considering that the atom, \vj) say, is to be removed from the atomic decomposition 
(TTTj) . and denoting the corresponding subspaces Vk\j and Wk\j, in order to span Wk\j the 
measurement vectors {|w^' 7 )}f =1 ^ - are modified according to the equation [2] 

i^-i^- '^y . < ■ j m ■ ■ /.. (17) 

Consequently, the coefficients in ffTTj) should be changed to 

c ^ = c "-^\¥> ' 1 /■•• (18) 

I If*) II 2 



5 Adaptive pursuit strategy for subspace selection 

Given a signal |/), we aim at finding the subspace Vk C V where one of the signal compo- 
nents lies. Let us stress once again that the problem arises from the impossibility of correctly 
computing the measurement vectors spanning the whole subspace W. Moreover, we have to 
face the fact that the corresponding signal component we want to represent is not available. 
What we know is that the available signal, |/), is expressible as the sum of two components 
\f) = l/v) + |/w x ) an d that there exists an unknown subspace Vk = s P an {|^)}|£i C V where 
{ii}^ is a set of K unknown indexes such that |/) = |/y) + |/w-l) = \fv K ) + l/w- 1 -)) with 

K 

i/v> = i/vk> = Ei^x ,w ^>- ( 19 ) 

i=i 

Hence, if the set of indexes {£i}fLi were given, one could construct the measurement vectors 
\wf-) in W K = sp&&{Pw\ v li}iLi an d obtain the component |/y) = \fv K ) fr° m P^]l - Unfortu- 
nately, in the problem we are addressing neither the set of indexes {li}f =l nor the component 
|/v) are given. Nevertheless, by applying Pw to both sides of (}T9]) we obtain 

K 

\M = \fw K ) = J2M( w ?W- ( 2 °) 

Denoting 1$ to the identity operator in S we have Pyy — Is — -Pw x • Thus, since the subspaces 
S and W ± are known, we do have access to the component |/w)- We can then look for the set 
of indexes to approximate this component as in (!20j) . 

Remark 3. Notice that the measures (wf\f) involved in (1191) and (1201) are the same. Therefore, 
by finding the representation of\f WK ) we have the information which is needed to obtain \ fv K ) 
from f|T9l) . Let us stress that the need to deal with \fw) o/so introduces the bad conditioned 
nature of the problem we are considering. Certainly, having access to the signal |/y) would 
imply that, provided that the spanning set {{vi)}^ were well conditioned, one could find the 
sparse approximation ffT^l) without difficulty. However, even when the conditioning of this 
spanning set is ideal (i.e. {Ifj)}*^ is an orthonormal basis for V) the fact that we need to deal 
with the projection \fw), which is sparse in the set {\ui) = Pw\vi)}iL\, introduces the difficulty 
we have to face. If the set {wj}*^ were well conditioned, the robust signal splitting could be 
obtained by a simple projection. The situation we are concerned with comprises the cases in 
which the whole set {-Uj}*^ is very bad conditioned but the projection ( |19p has a well conditioned 
representation in the subset {u^f =l C {wj}*^ we aim to find. 
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The proposed strategy for selecting the subset of atoms evolves by stepwise se- 

lection and is in line with the strategies in [23-25]. By fixing Pyv k , at iteration k + 1 we select 
the index i^+i such that ||JPyv|/) — Pw k+1 \f)\\ 2 is minimized. 

Proposition 1. Let us denote by J the set of indices {£\, . . . ,£k}- Given Wk = span^ue^}^, 
the index £k+i corresponding to the atom \ut k+1 ) for which \\Pyy\f) — Pw k+1 \f)\\ 2 is minimal is 
to be determined as 

4 + i = arg max ||| 7 „>||^0, (21) 

"GJ\J fc |||7n)|| 

with |7„) given in ( iTBl) . and Jk the set of indices that have been previously chosen to determine 
W fc . 

Proof. It readily follows since Av fc+1 |/> = AvJ/> + h Hf 1 and hence \\P w \f) - Av fc+1 |/> 1 1 2 = 
ll^w|/)|| 2 - ||AvJ/>|| 2 - K CT ! - Because P w \f) and P Wk \f) are fixed, ||P W |/) - Pw k+1 \f)\\ 2 
is minimized if tt^t ; ||7n|| 7^ is maximal over all n G J \ J^. □ 

The OBMP selection criterion [13] selects the index £k+i as the maximizer over n G J\ J& 

This condition was proposed in [13] based on the consistency principle [22,26]. Such a principle, 
introduced in [26] and extended in [22] , states that the reconstruction of a signal should be self- 
consistent in the sense that, if the approximation is measured with the same vectors, the same 
measures should be obtained. Accordingly, the above OBMP criterion was derived in [13] 
in order to select the measurement vector producing the maximum consistency error 

A = — E VkW ±f)\, with regard to a new measurement |w*ix)- However, since the 

measurement vectors are not normalized to unity, it is sensible to consider the consistency error 
relative to the corresponding vector norm 1 1 |iu^tj) 1 1, and select the index so as to maximize over 
k + 1 G J \Jk the relative consistency error 

llkOl ' IIK+l)ll#0 ' (22) 

In order to cancel this error, the new approximation is constructed accounting for the concomi- 
tant measurement vector. 

Property 3. The index £k+i satisfying (21\) maximizes over k + 1 G J \ Jk the relative consis- 
tency error (122]) 

Proof. Since for all vector |u;*+J) given in flB]) (w k k +l\E VhW ± = and = Wllk+i)^ 1 

we have 

A _ i(7 fc+ il/>l 

Hence, maximization of A over k + 1 G J \ Jk is equivalent to (12T1) . □ 

It is clear at this point that the forward selection of indices prescribed by proposition 
( 1211) is equivalent to selecting the indices by applying the subset selection criterion introduced 
in [23] (c.f. Theorem 1) on the projected signal Pw\f) using the dictionary In 
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the context of subspace selection for signal representation we have termed such a criterion 
Optimized Orthogonal Matching Pursuit (OOMP) [6]. 

The hypothesis that the computation of more than r measurement vectors becomes an 
ill posed problem enforces the OOMP selection of indices to stop if iteration r is reached. 
Nevertheless, the fact that the signal is assumed to be .fT-sparse, with K < r, does not imply 
that before (or at) iteration r one will always find the correct subspace. The r-value just 
indicates that it is not possible to continue with the forward selection, because the computations 
would become inaccurate and unstable. Hence, if the right solution was not yet found, one 
needs to implement a strategy accounting for the fact that it is not feasible to compute more 
than r measurement vectors. An adequate procedure is achieved by means of the swapping- 
based refinement to the OOMP approach introduced in [9]. As discussed below, it consists of 
interchanging already selected atoms with nonselected ones. 

Consider that at iteration r the correct subspace has not appeared yet and the selected 
indices are labeled by the r indices £\, . . . ,£ r . In order to choose the index of the atom that 
minimizes the norm of the residual error as passing from approximation PwJ/) to approxima- 
tion P-w r \j\f) we should fix the index of the atom to be deleted, £j say, as the one for which the 
quantity 

z = l,...,r. (23) 



is minimized [7, 9, 25, 27] . 

The process of eliminating one atom from the atomic decomposition ( TTTT) is called backward 
step while the process of adding one atom is called forward step. The forward selection criterion 
to choose the atom to replace the one eliminated in the previous step is accomplished by finding 
the index i = 1, . . . , r for which the functional 

e n = pHir , with \v n ) = \u n )-P Wr Ju n ), ||k>||^0 (24) 
1 1 r») 1 1 

is maximized. In our framework, using (flTl) . the projector Pw r \ 3 is computed as 



PW r \j — PWr 



Since Pvv r and \wj) are available, the computation of the sequence |z/ n ) in ( 124"1) is a simple 
operation. 

The swapping of pairs of atoms is repeated until the swapping operation, if carried out, would 
not decrease the approximation error. The implementation details for an effective realization 
of this process are given in [9], and MATLAB codes are available at [28]. Since there is no 
guarantee that at the end of the swapping of pairs of atoms the correct subspace has been 
found, the process can continue by increasing the number of atoms the swapping involves. At 
the second stage, in line with [29], we propose the swapping to be realized by the combinations 
of two backward steps followed by two forward steps, provided that the interchange of the two 
atoms improves the approximation error. If at the end of the second stage the right subspace 
has not yet been found, the number of atoms involved in the swapping is increased up to three 
and so on. Notice that if the number of atoms to be interchanged reaches the value r the whole 
process would repeat identically. This is avoided by initiating the new circle with a different 
initial atom. The above specified hypothesis ensures that the algorithm will stop when the 
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correct signal splitting has been found. At such a stage one has Pyj\f) = AvJ/) with W r 
spanned by the selected atoms £±, . . . , £ r . If the order K of sparseness of the signal is less than 
r a number of r — K coefficients in the atomic decomposition 

r r 
i=l i=l 

will have zero value. 

In the realistic case where the measurements are affected by errors, the proposed iterative 
process is to be stopped when the condition 

llAvl/>-Av r |/>ll<* 

is reached (where 5 should be determined by taking into account the errors of the data). 

6 Numerical simulations 
6.1 Impulsive noise filtering 

We extend here the example in [2] concerning the elimination of impulsive noise from the 
register of the motion of uncoupled damped harmonic oscillators. 

The n-th oscillator is characterized by a frequency of ^ Hertz and its motion, as a function 
of time, is given by the equation 

(t\x n ) = x n (t) = e~* cos(7rnt). (25) 
In [2] the register of system's motion was considered to be the signal 

. , „ . e cos [rent) r „ , , , 

W^>=g (i + o.7(n-ren - te|(u| (26) 

and the impulsive noise corrupting this signal was assumed to belongs to the subspace 

W ± = spanle" 100000 ^ - 0025 ^ 2 , j = 1, . . . ,400, t G [0, 1]}. (27) 

Due to the nature of the distribution of frequencies in (1261) . the corresponding oblique projector 
for filtering the impulsive noise from the register was recursively constructed by incrementing 
the frequency n one by one, until the projection of the noisy signal onto the span of the set 
{x n (t)}n =1 remained unaltered by increasing the number k of elements in the set up to some 
value. As remarked in [2], such a procedure is not always feasible. Let us consider for instance 
that the 100 oscillators have frequencies which are not restricted to the range [1, 100] but can 
be any integer number in the [1,405] interval. In this case the recursive construction of the 
projector by incrementing the frequency one by one is in general not possible, because the 
corresponding numerical calculations become ill posed before reaching the frequency n = 405. 
However, the oblique projector, along W 1 - given above, onto a subspace spanned by only 100 
functions x n (t), t G [0, 1] with n G [1,405], can be accurately calculated. Thus, by applying the 
techniques of the previous sections to find the right frequencies of the system, the cancellation 
of the impulsive noise in the register is rendered possible. This is illustrated by the numerical 
simulation described below. 
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The impulsive noise is simulated by randomly taking 200 pulses in (|2"Tj) . The system's 
motion is simulated as a linear combination of 100 functions ( 125]) the frequencies of which are 
taking, randomly, from the set [1,405]. The data are assumed to be known in single precision. 
The simulation was run 50 times and in all the cases the cancellation of the impulsive noise was 
successful. The left graph of Fig. 1 plots one of the realizations of the experiment (motion of 
the system plus impulsive noise vs time). The graph on the right depicts the result obtained by 
applying the proposed technique to the signal on the left (it coincides with the line representing 
the true signal). On the contrary, although the spaces fl27|) and span{x n (t), t £ [0, 1]}^=! are 
'theoretically' complementary, since the construction of the corresponding oblique projector is 
very ill posed, the projection fails to correctly separate the signals. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Figure 1: The left graph depicts the register of the motion of a system consisting of 100 damped 
harmonic oscillators, whose frequencies are integer numbers randomly taken from the interval 
[1,405], corrupted by 200 pulses randomly taken from the subspace W 1 - given in (J27J). The 
graph on the right depicts the result of filtering the signal of the previous graph by the strategy 
described in Section 5. The result coincides with the true motion of the simulated system. 



6.2 Application to the separation of a spectrum from blackbody 
radiation background 

We consider an hypothetical situation where the background is assumed to be produced by 
a linear combination of up to five blackbodies (e.g. stars) at temperatures of T\ = 3000K, 
T 2 = 3500K , T 3 = 4000K, T 4 = 4500K , and T 5 = 5000K. Hence, in this case the subspace W 1 
is defined as W L = span{?/j}f =1 , where i/i are functions of the wavelength A as given by 

C 

yi(X) = c± , Ci = 3.7419 x 10" 6 erg cm 2 s~\ C 2 = 1.4288cm K. 

A 5 (e4 - l) 

Any linear combination of these functions is an acceptable background. In our numerical 
experiment the background, y{\), is generated as y(X) = Y^=iViW- 

We simulate a spectrum, on a region of A ranging from zero to 3//m, by considering that it 
belongs to the cardinal cubic spline space on the interval [0, 3/xm] with separation b = 2 _4 /im 
between consecutive knots. Such a space can be spanned by a B-spline basis arising by translat- 
ing a prototype B-spline [30,31]. Different spectra are simulated by randomly drawing K = 70 
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functions {f^}™!, from the basis consisting of 483 functions, to generate the decomposition 
YliiLi c i v tn w ^h random coefficients q G [0, 1], i = 1 . . . , 70. A particular realization is de- 
picted in the bottom graph of Fig. 2. The simulated available signal is obtained by adding the 
background y{\) given above to the spectrum, and perturbing each data point with a normal 
distributed error of variance corresponding to a percentage of each data value. Let us remark 
that, although the background is not 'exactly' in the cardinal spline space of the spectrum, it 
has a very good representation in such a space. Hence, the calculation of an oblique projector 
onto the given spline space and along the space of the background is expected to be very badly 
conditioned. Indeed, for very small errors (variance of 10 _6 % of each data value) the separation 
of the spectrum from the background is not possible by an oblique projection onto the whole 
spline space. However, in a simulation of 100 different spectra, each of which consisting of 70 
randomly taken B-spline functions, the separation was successful in all the cases by applying 
the proposed approach. Thus, a more realistic situation was simulated by increasing the vari- 
ance of the error up to 1% of each data value. Also in this case the spectrum recovery was a 
complete success. In order to make evident the errors' effect, the variance was increased up 
to 5% of each data value. One of the simulations is plotted in the top graph of Fig. 3. The 
middle graph of the same figure shows both, the theoretical 'hidden' spectrum and the one 
recovered by the proposed approach. This graph is meant to illustrate a 'typical result', as in a 
run of the 100 simulations described above the reconstruction of the corresponding spectra was 
of similar quality. The visualization of the approximation quality is made clearer in the bottom 
graph of Fig. 3, where a portion of the previous graph (corresponding to the interval [0.5, 1]) is 
plotted. Here we can see that, as expected, the approximation (broken line) fails to reproduce 
the peaks of low intensity (the one at 0.65/im). This is of course understandable by comparing 
the intensity of the spectrum with the intensity of the data on the [0.5, 1] interval. In this 
region, a variance of 5% of intensity of each data point entails an uncertainty of more than one 
unit in the corresponding scale of intensity. Thus, one cannot expect to correctly spot peaks of 
intensity of the same order as the errors. On the other hand, some spurious peaks which are 
not in the true spectrum may also appear (note the small peak of negative intensity). However, 
on the whole we can confidently assert that the recovery of the simulated spectra is satisfactory 
even for significant error level. It is pertinent to point out that if one wished to avoid negative 
intensities one could penalize the selection of measurements leading to such negative values. In 
our framework this is implementable in a straightforward manner. For instance, in the forward 
selection procedure, at iteration k + 1 we select the index i^+i satisfying (1211 and it follows 
from (fTTI) and (fT5|) that the spectral intensity vector at this step is given as 

If \-lf \ ft U, \ (7fe+ll/) | i v (Tfe+iLQ / 98 n 
I/O - l/v ft ) - ^Wjjj^yjp + l"W> jjj^yjp (28) 

where |/y fc ) is the spectral intensity obtained in the previous iteration, E Vk the oblique projector 
onto the previously selected subspace, |/) the observed data vector and |7fc+i) (constructed as in 
f fTB"]) ) is to be determined in the selection process of the index £ k+ i according to the prescription 
of Proposition 1. Thus, restrictions on |/v fc+1 ) can be incorporated by disregarding the selected 
indices yielding unacceptable values of |/y fe+1 ). In the example we are discussing here, one can 
avoid negative values of intensity by disregarding those indices which satisfying (1211) do not 
fulfill the condition (A|/y fc ) = /y fc+1 (A) > for the values of A being considered. With the 
incorporation of this constrain the small negative pick in the middle graph of Fig. 3 disappears 
and the whole approximation improves. However, in some other realizations of the experiment 
when introducing the possibility constrains some small spurious picks of positive intensity 
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appear, yielding on the whole an approximation of quality comparable with the unconstrained 
one. The spurious small peaks (or the absence or small peaks present in the true spectrum) 
are consequence of the considerable uncertainty in the data. The fidelity of the approximation 
with the true spectrum (in all the realizations of the experiment) is improved by reducing the 
error of the data. 

7 Conclusion 

The construction of measurement vectors specially designed for separating signal components 
produced by phenomena of different nature was discussed. Assuming that the subspaces host- 
ing the signal components are given, the required measurement vectors should yield an oblique 
projection along one of the subspaces and onto the other. Considerations were restricted to 
those cases for which such subspaces are theoretically complementary, yet very close to each 
other, so that the construction of the measurement vectors for the whole space renders an ill 
posed problem. A recursive strategy for finding the right subspace to achieve the desired signal 
separation was then discussed. By recourse to numerical simulations it was illustrated that, 
provided that the signal is sparse in a spanning set of the signal subspace, the required signal 
splitting may be achieved by means of adaptive greedy techniques capable of searching for the 
required subspace while maintaining stability in the calculations. When tested in situations 
involving significant level of errors the proposed technique produced satisfactory results. There- 
fore, we are led to conclude that the framework for measurement design advanced here should 
be of assistance to a variety of applications where the discrimination of phenomena of different 
nature is required. 
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Figure 2: The top graph depicts theoretical data produced by the superposition of the spectrum 
shown in the bottom graph and black body radiation background. 
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Figure Caption 



Figure 3. The top graph has the same description as the top graph of Fig. 2, but the data 
are perturbed with zero mean normal distributed errors of variance corresponding to 5% of 
each data value. The continuous line in the middle graph represents the theoretical spectrum 
and the dotted line the spectrum obtained by the proposed approach from the data of the top 
graph. The bottom graph magnifies the region [0.5 1] in the previous graph. 
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